#ifndef BANDED_MATRIX_H
#define BANDED_MATRIX_H

#include<algorithm>
#include<functional>
#include<utility>
#include"layout.h"
#include"vector.h"
#include"tags.h"

namespace math{
template<class T>
class banded_matrix : public vector<T>
{
public:
    typedef vector<T> Array;
    typedef typename Array::difference_type difference_type;
    typedef typename Array::size_type size_type;
    typedef typename Array::value_type value_type;
    typedef typename Array::pointer pointer;
    typedef typename Array::const_pointer const_pointer;
    typedef typename Array::reference reference;
    typedef typename Array::const_reference const_reference;
    typedef typename Array::iterator iterator;
    typedef typename Array::const_iterator const_iterator;
    typedef typename Array::reverse_iterator reverse_iterator;
    typedef typename Array::const_reverse_iterator const_reverse_iterator;
    typedef matrix_tag category;

    template<class U>
    friend class banded_matrix;

    banded_matrix()
        :rows(1), columns(1), l(0), u(0),
		_rows(1), vector<T>(1)
    {}

	banded_matrix(size_type m, size_type n, size_type p, size_type q)
		:rows(m), columns(n), l(p), u(q), _rows(std::max(m, n)),
		vector<T>(std::max(m, n) * (p + q + 1))
	{
		assert(p < m);
		assert(q < n);
	}

	template<class U>
	banded_matrix(const banded_matrix<U>& m)
		:rows(m.rows), columns(m.columns),
		_rows(m._rows), vector<T>(m)
	{}

    inline size_type size1() const { return rows; }
    inline size_type size2() const { return columns; }
    inline size_type size_lower() const { return l; }
    inline size_type size_upper() const { return u; }

    inline T& operator()(size_type i, size_type j)
	{
		size_type ri = std::max(i, j);
		size_type ci = l + j - i;
		return data_[ci + ri * _rows];
	}

    inline const T& operator()(size_type i, size_type j) const
	{
		if((i > j && i - j <= l) || (i <= j && j - i <= u))
		{
			size_type ri = std::max(i, j);
			size_type ci = l + j - i;
			return data_[ci + ri * _rows];
		}
		else
			return zero_;
	}

private:
    size_type rows;
    size_type columns;
    size_type u;
    size_type l;
	size_type _rows;
    static value_type const zero_;
};

template<class T>
T const banded_matrix<T,Layout>::zero_ = T();

}


#endif
